      PROGRAM AUTOFLAT2                        
!
! Experimentalni verze 2.
! predchozi verze asi moc nefungovaly!!!
!
C Vyvoreni autoflat-fieldu z mnoha CCD snimku.
C
C FH, Rijen 1997
C FH, duben 1998
C FH, cerven 1998
! FH, srpen 1998
C
C Porad je to jen experimenatlni verze!
! Opatrnost nutna!
!
 
! Tato verze potrebuje krome obrazku taky nalezene hvezdy
! pro kazdy obrazek napriklad soubor .AP
    
       IMPLICIT NONE

       INTEGER MAXX,MAXY,NSIDE
       PARAMETER( MAXX = 1000, MAXY = 1000, NSIDE = 100 )

       CHARACTER*80 SB,BUF,SBap
       LOGICAL SIMPLE,EXTEND,ANYF,obloha(MAXX,MAXY),EXIST
       INTEGER I,J,L,JJ,eof,N,II
       INTEGER STATUS,BITPIX,NAXIS,NAXES(2) ! dimenze snimku
       INTEGER IMAGE(MAXX,MAXY),IBUF(MAXX*MAXY)
       INTEGER BLOCKSIZE,PCOUNT,GCOUNT
       INTEGER NAXISX,NAXESX(2),MEDIAN,VALUE
       integer nflat(maxx,maxy)
       REAL BSCALE,BZERO,FFLAT(MAXX,MAXY),x,y,r,smed,val,med
       
       EXTERNAL MEDIAN

c       DATA FFLAT/4 000 000*0.0/
       DATA NAXES/0,0/

       DO J = 1,MAXY
         DO I = 1,MAXX
           FFLAT(I,J) = 0.0
           nflat(i,j) = 0
         ENDDO
       ENDDO
       smed = 0.0

C otevreni souboru se jmeny souboru na vyrobu flat-fieldu
c       OPEN(1,FILE='af.in',STATUS='OLD')

       N = 0
! nacteni jmena souboru
   10  READ(*,'(A)',END=20) SB   
         if( sb .eq. ' ' ) goto 10 
         print *,sb(:50)
         n = n + 1
! otevreni souboru
         CALL FTOPEN(2,SB,0,BLOCKSIZE,STATUS)
         IF( STATUS.NE.0 )THEN
           PRINT *,'STATUS41=',STATUS
           STOP
         ENDIF    
! cteni parametru
         CALL FTGHPR(2,2,
     &        SIMPLE,BITPIX,NAXISX,NAXESX,PCOUNT,GCOUNT,EXTEND,STATUS)
         IF( STATUS.NE.0 )THEN
            PRINT *,'Spatne dimenze obrazku.'
            PRINT *,'STATUS48=',STATUS
            STOP
         ENDIF
         CALL FTGKYE(2,'BSCALE',BSCALE,BUF,STATUS)
         CALL FTGKYE(2,'BZERO',BZERO,BUF,STATUS)
! cteni obrazku
         L = 0
         DO J = 1,NAXESX(2)
           DO I = 1,NAXESX(1)
             L = L + 1
             CALL FTGPVJ(2,1,L,1,0,VALUE,ANYF,STATUS)
             IF( STATUS.NE. 0 .OR. ANYF )THEN
               PRINT *,'STATUS=',STATUS
c               STOP
             ENDIF             
             IMAGE(I,J) = VALUE
           ENDDO
         ENDDO

! vymazani hvezd z obrazku pomoci AP souboru:
         i = index(SB,'.fts')
         SBap = SB(:i-1)//'.AP'
         OPEN(15,file=sbap,status='old',iostat=i)
         if( i.ne.0 ) stop 'nenalezen AP'
         read(15,*)
         read(15,*)
         read(15,*)
         read(15,*)
         
! inicializace pomocneho pole
         DO j = 1,MAXy
           do i = 1,maxx
             obloha(i,j) = .true.
           enddo
         enddo
       
 110     call rdat(15,X,Y,eof)
!         print *,x,y
         if( eof .ne. 0 ) goto 120
         do i = -20,20
           do j = -20,20
             ii = i + x
             jj = j + y
             if( ii.ge.1 .and. ii.le.maxx 
     *      .and. jj.ge.1 .and. jj.le.maxy )then
             r = sqrt(real(i**2 + j**2))
             IF( r.le.15.0 ) obloha(ii,jj) = .false.
             endif
           enddo
         enddo
         goto 110
 120     close(15)


! z kazdeho 10-ho pixelu obrazku se udela lokalni median
        L = 0
        DO j = 1,Naxesx(2),5
          DO i = 1,Naxesx(1),5
            L = L + 1
            IBUF(L) = IMAGE(J,I)
          ENDDO
        ENDDO 
        MED = MEDIAN(L,IBUF,L/2)
        WRITE(*,*) '       Median=',MED
        smed = smed + med
!        MED = 1

      DO J = 1,NAXESX(2)
        DO I = 1,NAXESX(1)   
! hodnotu v obrazku pricte k vypocitavanemu flat-fieldu
           if( obloha(i,j) )then         
              val = IMAGE(I,J)
            FFLAT(I,J) = FFLAT(I,J) + val / MED
            nflat(i,j) = nflat(i,j) + 1
           endif
          
        ENDDO
      ENDDO
      CALL FTCLOS(2,STATUS)
      GOTO 10
   20 CONTINUE 

c       print *,i,j,fflat(i,j)
! konecny vypocet flat-fieldu
! konstatna 'smed' udava konecnou prumernou hodnotu (median)

       smed = smed / n
       DO J = 1,NAXESX(2)
         DO I = 1,NAXESX(1)
           if( nflat(i,j) .ne. 0 )then
            IMAGE(I,J) = smed*FFLAT(I,J)/nflat(i,j)
           else
            image(i,j) = smed
           endif
         ENDDO
       ENDDO

! cvicny vypocet medianu
        L = 0
        DO j = 1,Naxesx(2),1
          DO i = 1,Naxesx(1),1
            L = L + 1
            IBUF(L) = IMAGE(J,I)
          ENDDO
        ENDDO 
        MED = MEDIAN(L,IBUF,L/2)
        WRITE(*,*) '       Konecny median=',MED,' smed=',smed,n

       INQUIRE(FILE='autoflat.fts',EXIST=EXIST)
       IF( EXIST ) THEN
         OPEN(9,FILE='autoflat.fts')
         CLOSE(9,STATUS='DELETE')
       ENDIF
C Vytvoreni vystupniho obrazku
       CALL FTINIT(9,'autoflat.fts',1,STATUS)
       IF( STATUS.NE. 0 )THEN
         PRINT *,'STATUS103=',STATUS
         STOP
       ENDIF
C zapsani informace do hlavicky
       SIMPLE = .TRUE.
       BITPIX = 16
       NAXIS = 2
       NAXES(1) = NAXESX(1)
       NAXES(2) = NAXESX(2)
       CALL FTPHPS(9,BITPIX,NAXIS,NAXES,STATUS)
       CALL FTPKYS(9,'OBJECT','FLAT-FIELD','',STATUS)
       CALL FTPKYE(9,'BSCALE',BSCALE,10,'',STATUS)
       CALL FTPKYE(9,'BZERO',BZERO,10,'',STATUS)

C zapsani  dat
       CALL FTP2DJ(9,1,MAXX,NAXES(1),NAXES(2),IMAGE,STATUS)
       IF( STATUS.NE. 0 )THEN
         PRINT *,'STATUS120=',STATUS
         STOP
       ENDIF
C Zavreni souboru
       CALL FTCLOS(9,STATUS)

       print *,'Bye.'
       END       


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      INTEGER FUNCTION MEDIAN(N,A,K)
C find Kth element from array A
C According to procedure in
C N.Wirth: Algorithms + Data Structures = Programs, Prentice-Hall, 1975
C
      INTEGER N,K,L,R,I,J,W,X,A(*)

      L = 1
      R = N
10    IF( L.LT.R )THEN
        X = A(K)
        I = L
        J = R
30      IF( A(I).LT.X ) THEN
            I = I+1
            GOTO 30
          ENDIF
40        IF( X.LT.A(J) )THEN
            J = J-1
            GOTO 40
          ENDIF
          IF( I.LE.J )THEN
            W = A(I)
            A(I) = A(J)
            A(J) = W
            I = I+1
            J = J-1
          ENDIF
        IF( I.LE.J ) GOTO 30
        IF( J.LT.K ) L = I
        IF( K.LT.I ) R = J
        GOTO 10
      ENDIF
      MEDIAN = X
      END
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE RDAT(NUNIT,X,Y,EOF)
      
C  Cteni polohy ze souboru. Sam pozna, jestli jde 
C o typ AP nebo jiny a podle toho cte chyby jasnosti.

      INTEGER NUNIT,EOF
      REAL X,Y
      
      EOF=0

            READ(NUNIT,800,end=40) XX,X,Y
            READ(NUNIT,900,end=40) XX
            read(NUNIT,*,end=40)

800   FORMAT(F6.0,14F9.3) 
900   FORMAT(F13.0,2F6.2,F8.3,11F9.3)

      return

 40   EOF = -1
                             
      END


